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ABSTRACT 

This report summarizes the efforts to date in processing GPS measurements in High Earth Orbit (HEO) 
applications by the Colorado Center for Astrodynamks Research (CCAR). Two specific projects were 
conducted; initialization of the oibit propagation software, GEODE, using nominal orbital elements for the 
IMEX orbit, and processing of actual and simulated GPS data from the AMSAT satellite using a Doppler-only 
batch filter. 

CCAR has investigated a number of approaches for initialization of the GEODE oibit estimator with link a 
priori information. This document describes a batch solution approach that uses pseudorange or Doppler 
measurements collected over an orbital arc to compute an epoch state estimate The algorithm is based on 
limited orbital element knowledge from which a coarse estimate of satellite position and velocity can be 
determined and used to initialize GEODE. This algorithm assumes knowledge of nominal orbital elements, (a, e, 
i. OX fi) and uses a search on time of perigee passage (tp) to estimate the host satellite position within the orbit 
and the approximate receiver clock bias. Results of the method are shown for a simulation including la^s 
orbital uncertainties and measurement errors. 

In addition, CCAR has attempted to process GPS data from the AMSAT satellite to obtain an initial estimation 
of the orbit. Limited GPS data have been received lo date, with few' satellites tracked and no computed point 
solutions. Unknown variables in the received data have made computations of a precise orbit using the recovered 
pseudorange difficult. This document describes the Doppler-only batch approach used to compute the AMSAT 
orbit Both actual flight data from AMSAT. and simulated data generated using the Satellite Tool Kit and 
Goddard Space Flight Ct iter's Flight Simulator, were processed. Results ft each case and conclusion are 
presented. 
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INITIALIZATION ALGORITHM 
1. Introduction 

The GPS Enhanced Orbit Determination flight software (GEODE) was developed by Goddard 
Space Flight Center as a robust method for real-time mbit determination using the Global 
Positioning System. GEODE is a sequential Alter designed to process GPS pseudoranges and 
Doppler measurements from any number of satellites, with a dynamic model to determine mi 
estimate of the host vehicle orbit. The vehicle state is represented by position ?nd velocity 
coordinates in an earth-centered inertial frame (EC1). 

In low earth mbit (LEO) it is possible to initialize the GEODE state vector using a point 
solution i.e. an instantaneous solution for position and velocity based on four or more 
simultaneous pseudorange and Doppler measurements. If the state estimate is within a few 
kilometers of the true position, GEODE has been shown to consistently converge to an 
improved orbit solution within one orbit. This is not the case in high earth orbits or highly 
elliptical orbits, where there are rarely, if ever, sufficient satellites visible for a point solution. 
In these situations a different initialization approach is required. Furthermore, even in LEO, 
tite requirement for a point solution can delay the initialization of GEODE because of the need 
to perform a cold start/blind search for the satellite signals. 

The objective of our research has been to develop a robust method to initialize GEODE under 
the conditions expected in LEO and HEO. The basic assumption was that a batch approach 
was required to gather measurements over an arc of the orbit to produce an epoch state 
estimate sufficiently accurate to initialize GEODE. To provide a very robust approach, we 
assume only that the initialization algorithm has knowledge of the nominal orbital parameters 
for the vehicle. We then use typical launch vehicle oibit injection errors as the bounds for the 
errors in these starting points. We considered various parameterizations of the initialization 
batch and both search techniques and iterative batch schemes. 
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This document reviews the method used for this approach and presents the mathematical 
algorithms involved. Results from several test cases are also presented, as well as suggestions 
for future study. 

2. Approach 

The most obvious approach to initialization of GEODE would be to perform a batch solution 
for the position and velocity in ECI at a reference time. This has a straightforward 
relationship to the GPS pseudorange and Doppler measurements and directly produces the 
initial state estimate for GEODE. Unfortunately, there is no good way to construct this initial 
estimate of position and velocity. (Note: Earthbound users can assume an initial position at 
the center of the earth and achieve rapid convergence of a point solution approach). 
Furthermore, in order to tie together the time sequence of measurements for the batch, an orbit 
propagator must be used. The critical relationships among uncertainties in position and 
velocity elements within an orbit are not well captured by the ECI represet; tation. In fact, in 
our initial investigation we found that an initialization based upon position and velocity 
vectors is very sensitive to uncertainties and does not provide a robust approach. 

A better method, based upon an orbit dement representation, relies upon constraints of the 
orbital dynamics to narrow the range of possible initial conditions. This type of approach was 
selected based on the Tact tnat the nominal orbital elements arc well known and highly 
constrained by the launch trajectory. In addition, standard injection errors associated with 
these elements can be estimated based on the vehicle design and history. The advantage with 
this method is that the orbit errors provide far more geometric information than do the 
classical position and velocity error expressions. In particular, the angular position of the 
orbit (inclination, node, and argument of perigee), and orbital energy are very constrained by 
the launch; whereas the position of the vehicle within the orbit plane is not at all known. The 
values for the semi-major axis and eccentricity are somewhere in the middle, qualitatively, in 
terms of a priori knowledge. 

To describe the method clearly in the following section, we will use a classical orbit element 
representation and assume purely Keplerian motion. However, for the actual implementation 
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proposed for GEODE we outline the steps for an equinoctial representation using the dynamic 
propagation routines built into GEODE. 

2.1 Parameters 

For our development and initial testing of the initialization approach we parameterized the 
orbit using classical orbital elements (a, e, i, Q, to, and Tp). The first five elements define the 
geometry of the orbit with the nominal values for these parameters specified at launch. The 
final parameter, time of perigee passage (tp). together with the current time, prescribes the 
position of the vehicle within the orbit This value is not well known ahead of time. So, foe 
goal of foe initialization process is to identify the correct Tp and compute a position and 
velocity estimate for the satellite at foe filter start time. We considered both HEO and LEO 
orbits. 

2.2 Measurements 

We assume that the onboard GPS receiver acquires and tracks as many satellites as possible 
using a cold start or blind search technique. The receiver is assumed to form both pseudorange 
and Doppler measurements, and to collect the broadcast GPS satellite ephemeris data from the 
visible satellites. In HEO orbits, this visibility can be reduced to zero for extended periods of 
time as the host vehicle orbits above the GPS constellation. We assume that measurements 
are provided at 1 minute intervals, and consider solutions including pseudorange and/or 
Doppler. 

Time onboard the spacecraft is assumed to be known to within 1 second after acquisition and 
tracking of the first GPS satellite. A clock bias of up to 1 second may still be present, but the 
a priori estimate is set to zero. The stability of the clock is assumed to be 1 /10 10 . This means 
that over a 1 2 hour period the change in foe clock bias is within 5 ps. The overall clock bias 
is quite large (~ 3 x 1 0® m) but the drift error is less than 1 m/s. The large bias dominates the 
pseudorange r p duals, but we will show later that by comparing residuals for different 
satellites we can still compute a reliable initial condition using these measurements. If 
measurements from only one satellite are available it is not always possible to isolate the 
correct starting point in the orbit if there is a large clock bias. 
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23 Data Arc 

We considered various data arc lengths including 10 minutes, 100 minutes, and 600 minutes 
for the HEO orbit, and 10 minutes, 50 minutes, and 100 minutes for the LEO case, with 
measurements provided every 60 seconds. We also adjusted the starting point within the orbit 
to evaluate data arcs at perigee, apogee, and intermediate points. This provides a wide range 
of visibility and geometric conditions. 

2.4 Batch Solution 

The algorithm assumes knowledge of the orbital elements, (a, e, i, Q, to), and performs a 
search to estimate the remaining unknown - the location within the orbit, characterized by x p 
To process all the measurements in a batch solution, the nominal host vehicle elements are 
used to predict the vehicle position and velocity at each of the measurement times. For our 
initial investigation the propagation is though a purely Keplerian model with no perturbation 
effects. 

Combining the host vehicle position and velocity predictions with GPS satellite positions 
computed from the broadcast ephemerides, we compute the expected pseudorange and/or 
Dcppler measurements. These are compared to the measurements from the receiver and the 
residuals for the entire data arc are tallied. 

The characteristics of the measurement residuals for the batch indicate which t p is correct. 

For Doppler measurements and for PR measurements without clock biases, the RMS of the 
residuals is unambiguously smallest for the correct position within the orbit plane. In the 
presence of a large clock bias there is an offset in the residuals that prevents the use of a 
simple RMS evaluation. When observations are available from more than 1 satellite within 
the batch, the correct tp can be identified by the minimum standard deviation of the 
measurement residuals. This will be illustrated in the results shown in the following section. 
Thus, we compute the mean and standard deviation of the residuals for all measurements in 
the batch. The value of t p that minimizes the residual standard deviation locates the correct 
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host vehicle position within the orbit. The mean of the residuals for this ip provides a coarse 
estimate of the receiver clock bias. The position and velocity estimates can then be produced 
for the current time based upon the nominal elements and the optimal tp. 

The algorithm starts with a course search in x p , searching in increments of one tenth of an 
orbit. This search increment can then be refined to smaller resolutions to gain a better initial 
estimate of position. This is ultimately limited by the uncertainty in the nominal dements. To 
provide a general approach for different types of orbits, the current search implementation 
increments by fraction of an orbit period that correspond roughly to tens of minutes, one 
minute, and less than one second for both HEO and LEO orbits. In particular, we use 1/10, 
1/100, and 1/1 0000 of an mbit for the LEO case, 1/10, 1/1000 and 1/100000 of an orbit in foe 
HEO case. The figures in the following sections show that foe search space is well defined 
for data arcs of 100 min or longer eliminating foe possibility of searching in a false null 
region. For initialization near perigee, data arcs as short as 1 0 min are adequate. 

2.5 Computation of the Initial State Estimate 

Once the best Xpfx/) is found, the position and velocity of the host at a reference time are 
computed based on the assumed elements and x p \ The initial estimate of the clock bias may 
be taken as the mean of the pseudorange residuals for the minimum <Jsp residual solution. The 
initial covariance matrix may also be computed from the launch uncertainties. 

3. Simulation 

To evaluate the initialization approach prior to implementation in the GEODE code, we used a 
Matlab simulation environment. The components of the simulation are described in the 
following subsections. 

3.1 Truth Model 

We use M. Moreau's Matlab codes to establish the truth model for evaluating the algorithms. 
This includes 3 basic steps: 

a) Define host vehicle orbit (R,V, and orbit elements based on Keplerian model) and 
visibility conditions (antenna masking) 
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b) Define GPS satellite orbits (R,V computed from broadcast message file) and visibility 
conditions (antenna gain pattern) 

c) Compute range, range rate, and visibility (including antenna masking, gain patterns 
and Earth blockage) at 1 minute increments 

3.2 Measurement Model 

To produce measurements for the initialization algorithm we do the following: 

a) Select data batch (start and end times). Set the reference time to the start of the data 
batch. This insures that orbits are not propagated over very long intervals. 

b) Determine number of satellites observed and number of measurements per satellite for 
the entire arc. 

c) Produce error corrupted pseudorange and Doppler measurements. We add a large 
pseudorange bias (3e8 m) to all the measurements and a random error to each 
measurement (Gaussian with o=20m). Both the random error and bias are quite 
conservative. A random error is also added to each Doppler measurement (Gaussian 
with a = 5 m/s) 

3.3 Initialization algorithm 

a) Set up nominal orbit elements given standard injection errors. 

b) Search Tp over 1 orbit, incrementing in fractions of an orbit. 

For each x p : 

At each measurement time in the batch: 

i) . Compute position and velocity estimate based on nominal orbital elements 
and the current t p value. 

(a error, 0 error, *, flu* Xp , t meat) ^ (R» V) estimated at l™i 

ii) . Compute expected measurements and line of sights. Incorrect time 
estimates offset by one second are utilized for computing GPS satellite 
locations. 


Estimated LOS: 


LOS.. . = R, 


OPS truth 


-R 


SAT cjf 


Estimated Range: 

Pm =^T 
Estimated Range-rate: 



iii). Compute PR and/or Doppler residuals. These are accumulated over the 
batch for each * r 

Mcas Re sidual = [ p meas - p es i , Pmeas ~ Pest 1 
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SSR=SSR + MEAS 2 
P 

SSRdop = SSRdop + MEAS 2 

iv). The standard deviation and the mean of the measurement residuals are 
determined for each x p 

c) Compute R« t , V rat at reference time using best estimated x p and nominal orbital 
elements. 

3.4 Performance Analysis 

a. Evaluate metric for tp search to see if any false nulls occur 

b. Evaluate error in estimated R, V at the reference time to see if sufficient for 
GEODE initialization. 

3.5 Position and Velocity Calculation: 

Estimated orbital elements, (a, e, i, Cl, O), x p *), are converted to position and velocity estimates 
for direct comparison with the truth orbit. 


1 . Input orbital elements with injection errors, best estimate of x Pi reference time. 

2. Calculate (R, V)*,, at reference time: 

terror ' ^ error ' f p a , * Ci«u ) ^ (^/ > K )«/ 

3. Compare truth position and velocity with estimated position and velocity 

^ditt = Rirutk ~ Rest 

V - V - V 

rvz r lrHf h es! 

Compare truth position and velocity with estimated position and velocity in R1C 
coordinates 

a. Transform matrix from truth positions and velocities 
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^ truth r, r 

5 

truth t v; 


W- 


D yt' 

9X truth xys r tmth ire 

R y V 

lx truth tt r A m 


S=^x/? 


JVJ 


R(l) R(2) R(3) 
S(l) 5(2) 5(3) 
W(l) W(2) W(3) 


b. Compute position and velocity differences in RIC frame 


R *‘« ric 

V,„ 

"ft ric 


rirt 

* rtC 


^ rp .1 rs 
* nc 


/fxj 

r»i 


1 . Compute minimum position and velocity differences in XYZ, and RIC frames 
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4. Results 

To evaluate the algorithm we considered an ISS (LEO) and IMEX (HEO) orbit. Table 1 
summarizes the simulation parameters and tne errors considered. Two test scenarios with 
different orbit injection errors were considered for each type of orbit. We evaluated start 
times at perigee and apogee of the orbit using data arcs of 10 min, 100 min, and 600 min for 
the HEO orbit, and 10 min, 50 min and 100 min for the LEO orbit. 

Table 1. Simulation Parameters and Error Models 

Orbits 

ISS (LEO) (6415 km x 7091 km, 92 min period) 

IMEX (HEO) (6729 km x 42 164 km, 633 min period) 

Keplerian orbit propagation 
Orbit injection errors 
Scenario I 

perigee altitude: + 1 km 
apogee altitude: + 100 km 
inclination: +■ 1 degree 
nodes: + 1 degree 
Measurement errors 
Receiver clock bias - 0 s, 0. 1 s, 1 s 
Pseudorange - random Gaussian error (<j = 20) 

Doppler - random Gaussian error (a = 5) 

Data arcs 
HEO: 

Starting point - perigee (t = 1 or t = 633), apogee (t = 31 8 min), intermediate (t = 200 min) 

Arc length - 10 min, 100 min, 600 min 

LEO: 

Starting point- perigee (t = l), near apogee (t = 100), intermediate (t = 50) 

Arc length - 10 min, 50 min, 100 min LEO 


Scenario 2 

perigee altitude: + 1 km 
apogee altitude: + 10 km 
inclination: + 1/10 degree 
nodes: + 1/10 degree 
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Satellite visibility and altitude for the two orbits are shown in Figures 1 and 2, provided by 
M. Moreau’s visibility simulations. Data arcs and start times for each orbit type are also 
highlighted. 
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4.1 Simulation Resalts for ISS Initialization 

Figures 3 -5 show the T p search results for the ISS orbit given the larger injection errors 
presented in scenario I. Measurement residuals are computed for T p values spaceu at 0.1 orbit 
(approximately 10 minute) increments. The graphs each include curves for several starting 
points and data arc lengths. Figures 3 (a) and (b) show the r p dependence of the RMS 
residuals tor pseudorange and Doppler, respectively, in the case of no clock bias. Figure 4 
gives the RMS error of the pseudorange residuals with a I s clock error, and Figure 5 shows 
the mean and standard deviation of the pseudorange residuals for this case. 

From Figure 3 one can see that even data arcs as short as 10 min give a reliable minimum 
RMS residua! at the correct t p in the absence of a large clock bias. However, Figure 4 shows 
that the RMS residual is not reliable when large clock biases are present. Thus, we rely upon 
the standard deviation of the range residuals, shown in Figure 5b. The approximate clock bias 
value can be obtained as the mean value at t p * shown in Figure 5a. For example, the 100 min 
data arc starting at 50 gives a clock bias estimate of 0.0960 seconds (2.988 10e8 m). 

Figure 6 show s the RMS pseudorange and Doppler residuals for the narrower search on x p in 1 
s increments, over the region centered on the minimum value from the coarse search. The 
clock bias from the coarse search is included in the estimate. The minimum RMS value 
locates ip* to within 1 s. 

Using the a priori orbital elements and the derived T p \ the position and velocity of the host 
vehicle at the epoch are determined. Table 2 summarizes the errors in the epoch state estimate 
for »he ISS orbit given the injection errors of scenario l. Additionally, Table 3 summarizes the 
errors in the epoch state estimate from scenario 2 test runs. Values are given for x p resolution 
of I s up to 1 0 min 
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Figure 3a and 3b: RSS Residuals for 1SS Orbit with No Clock Bias. Graphs show 
multiple start times and data arcs for range residuals (top) and range rate residuals (bottom) as 
a function of the assumed t p . 
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Figure 4: RSS RttMub for ISS Orbit with 1 s Clock Bias. Graphs show multiple start 
times ami data arcs for range residuals as a function of the assumed tp. 
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Figure 5a and 5b; Mean and Standard Deviation of Range Residuals for ISS Orbit with 
1 s Clock Bias. Graphs show multiple start times and data arcs for range residuals as a 
(unction of the assumed t r Top is the mean of the residuals for each t P , and bottom graph is 
the standard deviation of the residuals for each T p . 
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Figure 6: RMS Range Residuals Over Narrow Search with Estimated Clock. Depicted 

are iterations on T p with intervals from 1/10 an orbit (approximately 10 minutes) to one 
second. 
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Table 2: ISS Position and velocity estimation errors for batch initialization. Initialization 
includes nominal orbit injection errors given in Table 1 for scenario 1 . Result r>~e shown for 
data arcs at perigee and apogee. Position and velocity errors tn radial, in-track, and cross- 
track directions are given for the end of the data arc ( 100 minutes). The correct time of 
perigee passage is 4342 $. 



Perigee 

1 ‘ 1 — — ■■ 

data arc (1038 measurements) 

Apogee 

data arc (962 measurements) 

Search 

Increment 

Tau 

Error (s) 

Position Error 
RIC (m) 

Velocity Error 
RIC (m/s) 

Tau 

Error* (s) 

Position Error 
RIC (m) 

Velocity Error 
RIC (m/s) 

0.1 orbit 
(552.3 s) 

126.32 

R = 3.76e4 
1 = 8.03e5 
C= 1.41e5 

R * 9.45e2 
1 = 87.37 
C = 70.91 

126.32 

R = 2.17e4 
1 = 1 28e6 
C = 1.3465 

R = 1.46e3 
1 » 1.57e2 
C = 64.138 

0 01 orbit 
(55.2 s) 

14.61 



R = 5.53e4 
1 = 334e4 
C = 1.48e5 

R = 14.17 
1 =35.13 
C = 51.25 

14.61 

R = 4 24e4 
1 = 4.16e5 
C= 1.409e5 

R =4.27 e2 
1 = 30.72 
C = 43.50 

0.0001 orbit 
(0.55 s) 


R = 5.84e4 
1 = 1.21eS 
C = f.49e5 

R = 72.96 
1 = 36.04 
C = 49.38 


10.149 

R = 4.26e4 
1 = 3.81e5 
C = 1.4le5 

R = 3.86e2 
1 = 28.29 
C = 42.64 


Table 3: ISS Position and velocity estimation errors for batch initialization. Initialization 
includes nominal (Kbit injection errors given in Table 1 for scenario 2. Position and velocity' 
errors in radial, in-track, and cross-track directions are given for the end of the data arc (10 
min). 



" — — * — "1 

Perigee data arc (96 measurements) 

Apogee data arc (1 13 measurements) 

Search 

Increment 

Tau 

Error (s) 

Position Error 
RIC (m) 

Velocity Error 
RIC (m/s) 

Tau 

Error * (s) 

Position Error 
RIC (m) 

Velocity Error 
RIC (m/s) 

0.1 orbit 
(552.3 s) 

82.068 

R = 4.08e4 
1 = 5.52e5 
C= 1.49e4 

R = 6.20e2 
1 = 2908 
C = 4.08 

82.068 

R = 2.99e3 
1 = 6.28e5 
C= 1.43e4 

R = 7.28e2 
1 = 3629 
C = 2.42 

0 01 orbit 
(55 2 s) 

l 

26.76 

R = 1.046e3 
1 = 2.39e5 
C = 1.51e4 

R = 1.59e2 
1 » 5.91 
C = 3.07 

26.76 

R = 8.72e3 
1 = 1.92e5 
C= i.45e4 

R = 2.29e2 
1 = 4.56 
C = 1 30 

0.00001 
orbit 
(0.55 s) 

8516 

R = 6.48e3 
1 = 2.57e3 
C= 1.52e4 

R = 7 42 
1 = 4.35 
C = 2.74 

3539 

R = 3.82e3 
I = 7.93e3 
C = 1 44e4 

R = 4.78 
l = 1.279 
C= 829 

i 
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Simulation Results for IMEX Initialization 

Similar results are given for the elliptical IMEX orbit. Figures 7 (a) and (b) show the * inge 
and Doppler RMS residuals for several start times and data arcs, without clock bias. t p values 
are spaced at 0. 1 orbit (63 min) increments. Figures 8 and 9 give the RSS, mean and standard 
deviation of the pscudorange residuals with a ! s clock error. 

Figures 7 (a) and (b) show that for a course search without clock bias, the correct x p * can be 
derived from the single minimum from both range and range-rate data. However, further 
examination at smaller search increments ( 1 0 minutes), the geometry of the range-rate data at 
apogee for data arcs of less than 600 minutes becomes insufficient for a single minimum 
search. 

As indicated in the ISS trials, the addition of a 1 second clock bias eliminates the possibility 
of using the RMS error to find the correct T p \ However, from Figure 9(b) one can see that the 
correct minimum can be located for long data arcs (approximately 600 minutes) and for start 
times near perigee. The geometry for short data arcs, especially near apogee, may still be 
insufficient for locating the correct Tp*. From Figure 9 (a), an approximate clock bias value of 
0.98505 sec (approximately 2.95 meters) can be determined. 

Using the nominal orbital elements and the derived t p * value, the position and velocity of the 
host vehicle at the epoch are determined and compared with truth values. Tables 4 and 5 
summarize the errors determined for the IMEX orbit in RIC coordinates for scenarios 1 and 2 
respectively. Uncertainties in the inclination and node elements negatively impact the 
accuracy of the x v * found within the search. This error in T P * significantly reduces the accuracy 
of the position and velocity estimates. 
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Figure 7a and 7b: Multiple Start time/Data Arc for IMEX case. Graphs show multiple 
start times and data arcs for range residuals (top) and range rate residuals (bottom) as a 
function of the assumed x p . 
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Figure 8: RSS Residuals for IMEX Orbit with 1 s Clock Bias. Graphs show multiple start 
times and data arcs for range residuals as a function of the assumed T p . 
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figure 9a: Mean of Range Residuals for IMEX Orbit with 1 s Clock Bias. Graph shows 
multiple start times and data arcs for range residuals as a function of the assumed t p . 
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Figure 9b: Standard Deviation of Range Residuals for 1MEX Orbit with 1 s Clock Bias. 
Graphs show multiple start times and data arcs for range residuals as a function of the 
assumed t p . 
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Table 4: Position and velocity estimation errors for batch initialization, initialization 
includes nominal orbit injection errors given in Table 1 for scenario 1. Results are shown for 
data arcs at perigee and apogee. Position and velocity errors in radial, in-track, and cross- 
track directions are given for the end of the data arc (100 minutes). The correct time of 
perigee passage is 38025 s. 



Perigee data arc (1038 measures -mts) 

Apogee data arc (982 measurements) 

Search 

Increment 

Tau 

Error (s) 

Position Error 
RIC (m) 

Velocity Error 
RIC (m/s) 

Tau 

Error* (s) 

Position Error 
RIC <m) 

Velocity Error 
RIC (m/s) 

0.1 orbit 
(3803.8 s) 

130.929 

R = 3.52e4 
1 = 8.33e5 
C = 4.99e5 

R = 72.38 
1 =82.76 
C = 60.9 

130.929 

R = 1.64e5 
1 = 1.13e6 
C = 5.9e5 

R = 31.51 
1 = 37.32 
C = 35.38 

0.01 orbit 
(380.38 s) 

16.46 

R - 2.54e5 
1 = 1.32e6 
C = 5.06e5 

R = 1 346e2 
1 = 80.49 
C = 59.80 

92.77 

R = 1.2465 
1 = 1.19e6 
C = 5.95e5 

R = 41.33 
1 = 37.61 
C = 35.63 

0.00001 

orbit 

(0.380 s) 

13.36 

R = 2.61 e5 
1 = 1.36eg 
C = 5.06e4 

R = 1.36e2 
1 = 80.43 
C = 59.77 

85.14 

R = 1.15e5 
I = 1.20e6 
C = 5.93e5 

R = 43.30 
! = 37.67 
C = 35.76 


Table 5: Position and velocity estimation errors for batch initialization. Initialization 
includes nominal orbit injection errcrs of scenario 2, found in Table 1 . 



Perigee data arc (1038 measurements) 

Apogee data arc (982 measurements) 

Search 

Increment 

Tau 

Error (s) 

Position Error 
RIC <m) 


Tau 

Error* (s) 

Position Error 
RIC (m) 

Velocity Error 
RIC (m/s) 

01 orbit 
(3803.8 s) 

25.837 

R = 3.326e4 
I = 5.l97e4 
C * 4.956e4 

R = .457 
1 = 8.41 
C = 6.12 

25.8 

R = 3.10e4 
1 = 9.14e4 
C = 5.99e4 

R = 80 03 
1 = 93.61 
C = 3.53 

0.01 orbit 
(360.38 s) 

25.8 

R = 3.31e4 
1 = 5.20e4 
C = 4.95e4 

R = .466 
1 = 8.40 
C = 6.10 

1221 

R = 8.01 e3 
1 = 1.1 5e5 
C = 5.98e4 

R = 9.97 
1 = 3.64 
C = 3.55 

0.00001 
crbit 
(0.380 s) 

1594 

R - 6.93e3 
1 = 7.73e4 
C = 4.96e4 

R = 5,97 
1 =8.39 
C = 6.11 

8.714 

R = 1 30e4 
1 = 1.21e5 
C = 5,89e4 


R = 4.48 
1 = 3.62 
C = 3.54 


2 b 
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Initialization Results Analysis 

The initialization search procedure described within this memo was shown to be effective in 
location the position of a satellite within its orbit, based on GPS observations. The accuracy of 
the location is limited by the accuracy of the assumed orbit. Thus, if the errors in the nominal 
orbit elements are sufficiently small to initialized GEODE, this procedure allows you to solve 
for locations within the orbit to the same level of accuracy. As the estimates of the other 
elements have not been adjusted, the position and velocity estimates cannot be better than the 
a priori. For example, if the error in the assumed perigee is 10 km, the resulting perigee 
estimate cannot be better than 1 0 km, but we can determine if the satellite is at perigee. 

Subsequently, based on these results, an additional step is required to develop a procedure to 
adjust all the element estimates based on the data so that the final position and velocity 
estimates are with in the convergence regime of GEODE. 
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AMSAT DOPPLER PROCESSING 


I. Introduction 

The AMSAT OSCAR-40 (AO-40) spacecraft was launched into a highly dliptical cibit 
(H tO) in 2001 . Although primarily designed to support amateur radio experiments, secondary 
experiments have been conducted to determine the feasibility of using GPS measurements for 
position and attitude determination in HEO. 1 he satellite is in a low inclination, 1000 x 
58,800 km altitude orbit, with an orbital period of approximately 19 hows ACMO uses an 
existing LEO GPS receiver, the Trimble TANS Vector, which poses sipuficant operating 
limitations at high altitudes above the GPS constellation. However, the AO-40 experiment has 
already proven successful, achieving autonomous acquisition and tracking of GPS signals 
throughout the orbit including signals with high dynamics around perigee and weak signals 
near apogee 

The TaNS Vector reports sev eral GPS observed variables, including code phase, carrier 
phase, and Doppler. Although the receiver is capable of returning position, velocity, and dock 
solutions when four or more GPS satellites are simultaneously tracked, no point solutions 
have been received to date. The primary tracking observable, pseudorange, can be constructed 
from the repe'.ed code phase. However, the lack of point solutions creates several unknowns, 
including satellite positions and clock biases, which make pseudorange recovery difficult. 

The objective of our research was to utilize the returned Doppler measurements for post- 
processing orbit determination. Although Doppler solutions are generally less accurate, the 
Doppler measurements are less sensitive to the effects of the receiver, and may provide an 
iniual estimation of the orbit from which to initialize other high fidelity post-processing 
schemes such as GEODE. Actual flight data from the receiver onboard AMSAT were 
processed in a Doppler-only batch filter. In addition, to further evaluate the performance of 
the filter and characterize unknowns within the flight data, simulated data were also processed 
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in the i.tter. This document reviews die Doppler-onlv batch algorithm, and presents the 
results bom actual and simulated data processmr 


2. Batch Approach 

We implemented a traditional hatch filter to process the Doppler measurements, the basis of 
which is similar to the batch described for the initialization algorithm The filter is utilized to 
process AMSAT flight data, consisting of returned Doppler data, as well as simulated data, 
based on generated satellite positions and velocities 

This approach assumes an initial ’truth’ vehicle state, and propagates through the orbit to each 
of the measurement times. Initially, this propagation uses a purely Keplenan model without 
considering perturbation effects. Estimated range-rates are computed from the propagated 
satellite positions and velocities, and the GPS satellite positions determined from the 
broadcast ephemeredes. The estimated range-rate is compared with either the measured 
Doppler converted to range-rate, or a simulated range-rate to produce measurement residuals 
for the entire arc of data. A least squares solution is performed to produce an updated estimate 
of initial satellite position and velocity, and iterations are performed until the upd ited estimate 
and the initial truth match withm tolerances. This simple filter was additionally updated to 
consider measurement w eighting based on a priori information. The results of both the simple 
filter and the weighted batch filter are presented and compared. 

2.1 State Parameters 

The vehicle state is comprised of position, v elocity, and frequency bias. 

X - [/?, R, R : R, R, R ; f 
L 

Initially, time tag error. TE. was also included in the estimated state. However, further 
analysis of the Doppler measurement partials indicated that the sensitivity' of range-rate to 
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l 
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time tag error is low ( 1 sec TE - 10 cm/s range-rate). With time tag errors expected to be 
fairly small, adequate solutions should be achievable without TE estimation, and therefore it 
was not considered in the final solution. 

2.2 Measurements 

GPS Receiver 1 onboard AO-40 was first operated from September 25* through November 
2 nd , 2001. This receiver uses a blind search technique, sequentially searching though the 
possible GPS satellite PRNs to acquire and track GPS satellites. The most promising batch of 
data was obtained on October 5 th , cm a pass near perigee were 4 satellites appear to be tracked 
simultaneously, although point solutions were not returned. The receiver forms Doppler 
measurements for each satellite tracked, which are provided at approximately half-second 
interv als. A 1 5 minute arc of data with 5 total PRNs tracked around the perigee pass was 
considered. The observed Doppler measurements were converted to range-rate measurements 
for processing in the filter. 

Simulated observations were also constructed for the same arc of data using satellite positions 
and velocities generated in STK. The measured range-rate was constructed from a non- 
linearized true range-rate model plus a random error. 


3. Data Processing 

The batch algoritfun used a MATLAB simulation environment to post-process both the flight 
data and the simulated STK data, the components of which are described below. 

3 . 1 Reference T ra jectory 

The reference trajectory for AO-40 was generated from a NORAD Two-1 ine Ephemeris 
(TLE), tor October 5 ;h , 2001 using STK’s MSGP4 (Merged Simplified General 
Perturbations) propagator. This high fidelity propagator considers secular and periodic 
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variations due to J2, solar and lunar effects, and atmospheric drag. The position and velocity 
at the initial measurement time are pulled from this reference trajectory and used as the initial 
state estimate in the batch propagator. 

3.2 Measurement models 


The batch filter utilizes the same process to generate computed measurements for both 
AMSAT flight data and simulated data scenarios. However, computations for observed 
measurements for each scenario differ slightly, and will be discussed separately. 


3.2.1 Computed Measurements 

a) Initial state estimate from the STK reference trajectory is propagated to each of the 
measurement times using a 2-body propagator 

b) The estimated positions and velocities, and GPS positions and velocities are used to 
generate an estimated range-rate at each of the measurement times, based on a non- 
linear Doppler model. 

cj Estimated range-rates are calculated similarly for both flight data and simulated data 
scenarios 

3.2.2 Observed Measurements 

AMSAT Flight Data 

a) The batch filter selects the Doppler measurement and converts to a range-rate 
measurement by multiplying by the negative LI wavelength. 

STK Simulated Data 

b) The measured range-rate is derived from the reference trajectory. 

c) The reference positions and velocities at each measurement time are used with the 
GPS positions and velocities to generate a ‘true range-rate’. 

dl A random noise (few m/s) is then added to each true measurement to create a 
measured range-rate. 

3.3 Doppler-only Batch Algorithm 

a) Determine reference trajectory in using a TLE in STK 

b) Compute AMSAT R/V at the measurement time based on the a priori R, v V„ from 
reference trajectory propagated in the 2-body propagation algorithm 

c) Calculate the estimated range-rates 
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Estimated LOS: 

LOSe, = Rcrs- Ramsat 


Estimated Range: 
Estimated Range-rate 



d) Determine measured range-rates 
a. For AMSAT flight data: 

Measured range-rate: 


b. For simulated data 

Trje range-rate: 



Measured range-rate: 

= Prw + random noise 


c) Compute Doppler-only measurement residuals 


0 Form H matrix 


dp_ , 

dR dr 


-V 


AMSAT 


♦ b • 

r eft 


Los al 

Pet > 


Pr, 


ext 

Pn, 


l 


J 
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g) Accumulate measurement residuals and H matrices for all measurement times. 

h) Solve for the state correction, dx. 

dy - H*dx 

i) Update the initial position velocity with dx , iterate. 


4. Data Processiag and Results 

AMSAT FLIGHT DATA 


Actual Doppler data consisting of a fifteen minute arc near perigee were processed in the 
batch filter, summarized in Table 1 . 


Table 1: Flight Data 


Date 

Start SOW 

End SOW 

Data Rate 

Satellites 

Notes 

10.05.01 

493194.2513 



494096.75 



- 0.5 seconds 

3, 11,22, 25,31 

Near perigee, 

4 simultaneous 


An initial state estimate from STK was propagated to the measurement times using the 
Kepi cn an propagator. Other perturbations were not considered. Additionally, a priori 

weighting of the data was not considered. The batch filter estimated values of , and 

the observed minus computed measurement residuals were examined over several iterations to 
evaluate the performance of the filter. 

Initially, examination of the residuals for the first two iterations indicates that the filter is 
adjusting the estimates properly . The residuals begin to converge to a zero value as expected, 
shown in Figure I . In addition, the filter appears to attempt to converge over the next 10 

iterations. The values of appear to converge to appropriate values, while the state 

correction value, dx , for each variable converges on zero as anticipated, shown in Figures 2. 3 
and 4 respectively. Figure 5 illustrates the residuals over these 10 iterations. 
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Although evaluations of the state parameters indicate that the filter is slowly trending toward 
convergence, examination of the residuals reveal potential problems. Although the residuals 
appear to center around a zero value, the filter appears unable to converge the measurements 
from each separate GPS satellite. Further, additional iterations do not continue the trend 
towards convergence, and a suitable solution is not obtained. Instead the filter appears to 
become unstable, with the residuals becoming very large at approximately 20 iterations, as 
shown in Figure 6. 


38 




CCAR 


Final Report 



Figure 6: Measurement Residuals {m/s), 21 iterations. 


Initially, we considered the possibility that improper measurement weighting contributed to 
the lack of convergence and ultimate failure of the filter in the previous scenario. To examine 
this further, the filter was expanded to incorporate a priori covariance and weighting 
estimations. The general solution of this weighted least-squares filter is shown below. 

Solve L * v o = N for dx 

’• vhere L = Pf' + ^ HWH and N - Pf { *dx + y 

P - covariance matrix 
W - weighting matrix 
x„- initial state deviation, or correction 

However, multiple attempts to adjust the weighting and a priori covariance did not appear to 
correct the filter divergence, as shown in Figure 7. 
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Figure 7: Measurement Residuals for Weighted Least-Squares Filter (m/s), 20 iterations 


Several unknowns exist with the flight data that could contribute to the failure of the filter. 
These include bad data related to one or more PRNs, poor geometry of the GPS satellites, not 
enough data in the batch, and unmodeled noise within the data. To examine these possibilities, 
we constructed simulated data with fewer unknown error sources for the same batch. The 
results of processing these data are discussed below. 


STK SIMULATED DATA 


Using the positions and velocities from the reference trajectory, truth range-rate measurements 
were computed. To these we added a measurement error of less than la to generate the 
‘observed’ measurements. An offset of 1% in semi-major axis and eccentricity was added to 
the initial state estimate from the reference trajectory from which the estimated AMSAT 
positions, velocities and range-rate measurements were calculated. Additional parameters. 
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such as number of satellites, data rate and total batch size remained the same as in the flight 
data scenario. 

Figure 8 illustrates the measurement residuals for the simulated data after 20 iterations. 
Figures 9(a) and (b) show the state corrections, dR and dV, e or each iteration. It is clear that 
the filter is abie to process the data successfully, converging on a suitable solution after 
approximately eight iterations. 



Figure 8: Measurement Residuals for Simulated Data (sn/s), 20 iterations 
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Table 2 contains the initial offset values, the final values for the state estimate at epoch, R<*V 0 , 
as well as the difference between the final solution and truth alter 20 iterations. 


Table 2: Siaalated Data Result* 


IffffWBTTWl 

■EEETSn 

Estimated 


Trath-Est 

Trath-Est 

(1% a, 1% e) 

(1% a, 1% e) 

laitial 

NESS 1 

Position 

Velocity 

Position (m) 

Velocity (m/s) 

Position (R,) 

ISHEzBi 

(m) 

(m/s) 

-1.694 e5 

-4.4 e2 

-5.807e6 

6.048e3 

34.76 

.491 

-1.35 c5 

5.6 e2 

-4.607e6 

-7.67c3 

86.60 

.831 

- 093 e5 

-0.2 e2 

-0.32 Ie6 

0.946c3 

-397.3 

-.257 


These results clarify several concerns First, it verifies the filters ability to process data 
correctly, and rules out any errors in software or the algorithm More importantly, it indicates 
that both the geometry of the GPS satellites, as well as the amount of data processed is 
sufficient for filter stability and convergence. Given this, it seems likely that the filter fails on 
the actual flight data because of unknown and therefore poorly modeled errors or noise in the 
Doppler data. It is interesting to note that the filter will fail with the simulated data if too 
much measurement noise (more than a few m s) is added. Actual flight data are expected to 
have more measurement noise, therefore we reasonably believe this is one specific cause of 
filter failure on the AMSAT data. 


5. Fata re Studies 

Goddard Space Flight Center has generated additional data using an actual TANS Vector 
receiver in a high fidelity flight simulator for the Landsat 7 orbit (LEO). Again, being 
simulated, these data have fewer unknowns within the scenario, and therefore can be used to 
examine the function of the receiver itself and the quality of the data generated. Preliminary 
attempts to process these data tn the Doppler-only batch were performed to further determine 
what clemerts within the flight data from this receiver create divergence in the filter. 

Simple studies using STK indicate that effects from the oblatencss of the Earth, J2, are more 
pronounced with a LEO orbit than with a HEO, as expected. Figure 10 compares the Landsat 
7 orbit propagated with a 2-body propagation verses propagated with J2 effects included. The 
errors resulting from a 2 -body propagation in LEO indicate that the Keplerian propagator 
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measurements, the in-track position of the estimated orbit steps and then rapidly diverges 
from the truth m-track position. The same trend cart he seen in the radial velocity 
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Figure 1 1 : Truth Orbit - Estimated Orbit (RIC), LS7 fm, tali) 


This propagator has successfully generated orbits for other similar orbit scenarios, so it is 
unclear why these trends are seen with the Landsat 7 orbit. It is believed that there may be a 
coordinate system difference between truth and the estimated orbit but this has not been 
confirmed. Successful testing with this data set could be extremely useful in understanding 
the data returned from the T ANS Vector, and for further post-processing attempts of data 
from AMSAT. 


Results Analysis 

The Doppler-only batch procedure described was shown to be effective in post-processing 
orbit determination for the simulated AMSAT orbit, based on generated measurements and 
visible GPS satellite information. AMSAT flight data have yet to be filtered successfully, 
believed in part to be due to noise and additional unknown errors in the flight environment 
and with the TANS Vector receiver itself Successful testing with the Landsat 7 data set could 
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he extremely useful in understanding the data returned from the T ANS Vector, and in further 
post-processing attempts of data from AMSAT. 
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